// ==================================================================================
// Shared Genomics Project MPI Codebase
// Version 1.0 30/04/2010
//
// (c) 2010 University of Manchester all rights reserved
//
// This file is distributed under the GNU General Public License, Version 2.  
// Please see the file COPYING.txt for more details
// ==================================================================================

/*!
\file 
\brief Implementation of the PMATHS Library 
*/

#include <ctype.h>
#include <float.h>
#include <string.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#include "cdflib.h"
#include "pmaths.h"
#include "utils.h"
#include "mt19937ar.h"

#define MAX(a, b) (((a) > (b)) ? (a) : (b))
#define MIN(a, b) (((a) < (b)) ? (a) : (b))
#define SQR(a) a * a
#define SIGN(a, b) b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a)

#ifndef MESSAGE_BUFFER_LENGTH
#define MESSAGE_BUFFER_LENGTH 1024
#endif

#define FPMIN 1.0e-30
#define vif_threshold  (50.00)

BOOL realnum(double d) {
	double zero = 0;
	if (d != d || d == 1/zero || d == -1/zero) 
		return FALSE;
	else return TRUE;
} 

double chiprobP(double x, double df) {
	double p, q;
	int st = 0; // error variable
	int w = 1; // function variable
	double bnd = 1; // boundary function

	if (!realnum(x)) return -9;

	// NCP is set to 0
	cdfchi(&w,&p,&q,&x,&df,&st,&bnd);

	// Check status
	if (st != 0 ) return -9;

	// Return p-value
	return q;
} 

int permute(int *a, int nA)
{
	int i, j;
	if (a == NULL || nA == 0) return 0;

	// generate random permutation of 0..n-1
	// where nA is the size of the input array.
	for(i = 0; i < nA; i++ )
		a[i] = i;

	for(j = 1; j < nA; j++ )
	{
		long int pos = get_rand(j+1);
		long int tmp = a[j];

		a[j] = a[pos];
		a[pos] = tmp;
	}

	return 1;
}

int get_rand(int n) {
	int r;

	r = (int) (n * genrand_real2());
	if (r == n) r--;

	return r;
}

double gammln(double xx)
{
	double x, y, tmp, ser;
	static double cof[6]={ 76.18009172947146, -86.50532032941677,
		24.01409824083091, -1.231739572450155,
		0.1208650973866179e-2, -0.5395239384953e-5};
	int j;
	static double five_five = 5.5;

	y=x=xx;
	tmp=  x +  five_five;
	tmp -= (x+ 0.5)* log(tmp);
	ser= 1.000000000190015;
	for (j=0; j<=5; j++) ser += cof[j]/++y;
	return -tmp+  log(2.5066282746310005*ser/x);
}

double betacf(const double a, const double b, const double x)
{
	const int MAXIT = 100;
	const double EPS = 3e-7;

	int m,m2;
	double aa,c,d,del,h,qab,qam,qap;

	qab=a+b;
	qap=a+  1.0;
	qam=a-  1.0;
	c= 1.0;
	d= 1.0-qab*x/qap;
	if (fabs(d) < FPMIN) d=FPMIN;
	d=1.0/d;
	h=d;
	for (m=1;m<=MAXIT;m++) {
		m2=2*m;
		aa=m*(b-m)*x/((qam+m2)*(a+m2));
		d=1.0+aa*d;
		if (fabs(d) < FPMIN) d=FPMIN;
		c=1.0+aa/c;
		if (fabs(c) < FPMIN) c=FPMIN;
		d=1.0/d;
		h *= d*c;
		aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
		d=1.0+aa*d;
		if (fabs(d) < FPMIN) d=FPMIN;
		c=1.0+aa/c;
		if (fabs(c) < FPMIN) c=FPMIN;
		d=1.0/d;
		del=d*c;
		h *= del;
		if (fabs(del-1.0) <= EPS) break;
	}
	if (m > MAXIT) return 0.0;
	return h;
}

double betai(const double a, const double b, const double x)
{
	double bt;

	if (x < 0.0 || x > 1.0) return 0;
	if (x == 0.0 || x == 1.0) bt=0.0;
	else
		bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x));
	if (x < (a+1.0)/(a+b+2.0))
		return bt*betacf(a,b,x)/a;
	else
		return 1.0-bt*betacf(b,a,1.0-x)/b;
} 

double pT(double T, double df)
{
	// Old NR routine
	double p = betai(0.5*df,0.5,df/(df+T*T)); 

	if ( p <= 1e-20 ) p = 1e-20;
	else if (p>=1) p=1; 

	return p;
} 

// Inverse normal distribution

/*
* Lower tail quantile for standard normal distribution function.
*
* This function returns an approximation of the inverse cumulative
* standard normal distribution function.  I.e., given P, it returns
* an approximation to the X satisfying P = Pr{Z <= X} where Z is a
* random variable from the standard normal distribution.
*
* The algorithm uses a minimax approximation by rational functions
* and the result has a relative error whose absolute value is less
* than 1.15e-9.
*
* Author:      Peter J. Acklam
* Time-stamp:  2002-06-09 18:45:44 +0200
* E-mail:      jacklam@math.uio.no
* WWW URL:     http://www.math.uio.no/~jacklam
*
* C implementation adapted from Peter's Perl version
*/

// Coefficients in rational approximations. 
static const double _a[] =
{
	-3.969683028665376e+01,
	2.209460984245205e+02,
	-2.759285104469687e+02,
	1.383577518672690e+02,
	-3.066479806614716e+01,
	2.506628277459239e+00
};

static const double _b[] =
{
	-5.447609879822406e+01,
	1.615858368580409e+02,
	-1.556989798598866e+02,
	6.680131188771972e+01,
	-1.328068155288572e+01
};

static const double _c[] =
{
	-7.784894002430293e-03,
	-3.223964580411365e-01,
	-2.400758277161838e+00,
	-2.549732539343734e+00,
	4.374664141464968e+00,
	2.938163982698783e+00
};

static const double _d[] =
{
	7.784695709041462e-03,
	3.224671290700398e-01,
	2.445134137142996e+00,
	3.754408661907416e+00
};

#define LOW 0.02425
#define HIGH 0.97575

double ltqnorm(double p)
{
	double q, r;

	errno = 0;

	if (p < 0 || p > 1)
	{
		errno = EDOM;
		return 0.0;
	}
	else if (p == 0)
	{
		errno = ERANGE;
		return -HUGE_VAL; // minus "infinity";
	}
	else if (p == 1)
	{
		errno = ERANGE;
		return HUGE_VAL; // "infinity";
	}
	else if (p < LOW)
	{
		// Rational approximation for lower region
		q = sqrt(-2*log(p));
		return (((((_c[0]*q+_c[1])*q+_c[2])*q+_c[3])*q+_c[4])*q+_c[5]) /
			((((_d[0]*q+_d[1])*q+_d[2])*q+_d[3])*q+1);
	}
	else if (p > HIGH)
	{
		// Rational approximation for upper region
		q  = sqrt(-2*log(1-p));
		return -(((((_c[0]*q+_c[1])*q+_c[2])*q+_c[3])*q+_c[4])*q+_c[5]) /
			((((_d[0]*q+_d[1])*q+_d[2])*q+_d[3])*q+1);
	}
	else
	{
		// Rational approximation for central region
		q = p - 0.5;
		r = q*q;
		return (((((_a[0]*r+_a[1])*r+_a[2])*r+_a[3])*r+_a[4])*r+_a[5])*q /
			(((((_b[0]*r+_b[1])*r+_b[2])*r+_b[3])*r+_b[4])*r+1);
	}
}

double pythag(const double a, const double b)
{
	double absa,absb;

	absa=fabs(a);
	absb=fabs(b);
	if (absa > absb) return absa* sqrt(1.0+SQR(absb/absa));
	else return (absb ==  0.0 ?  0.0 : absb*sqrt(1.0+SQR(absa/absb)));
}

/*
// This function implements an exact SNP test of Hardy-Weinberg
// Equilibrium as described in Wigginton, JE, Cutler, DJ, and
// Abecasis, GR (2005) A Note on Exact Tests of Hardy-Weinberg
// Equilibrium. American Journal of Human Genetics. 76: 000 - 000
//
// Written by Jan Wigginton
*/
double SNPHWE(int obs_hets, int obs_hom1, int obs_hom2)
{
	int obs_homc, obs_homr, rare_copies, genotypes, i, mid;
	VECTOR_T het_probs = NULL; 
	double sum, p_hwe = 0.0;
	int curr_hets, curr_homr, curr_homc;

	if (obs_hom1 + obs_hom2 + obs_hets == 0 ) return 1;

	// "Internal error: negative count in HWE test.
	if (obs_hom1 < 0 || obs_hom2 < 0 || obs_hets < 0) 
		return -1;		
	
	obs_homc = obs_hom1 < obs_hom2 ? obs_hom2 : obs_hom1;
	obs_homr = obs_hom1 < obs_hom2 ? obs_hom1 : obs_hom2;

	rare_copies = 2 * obs_homr + obs_hets;
	genotypes   = obs_hets + obs_homc + obs_homr;

	het_probs = create_vector(rare_copies + 1, 0);
	if (het_probs == NULL) return -1;

	// start at midpoint 
	mid = rare_copies * (2 * genotypes - rare_copies) / (2 * genotypes);

	// check to ensure that midpoint and rare alleles have same parity
	if ((rare_copies & 1) ^ (mid & 1))
		mid++;

	curr_hets = mid;
	curr_homr = (rare_copies - mid) / 2;
	curr_homc = genotypes - curr_hets - curr_homr;

	het_probs[mid] = 1.0;
	sum = het_probs[mid];
	for (curr_hets = mid; curr_hets > 1; curr_hets -= 2) {
		het_probs[curr_hets - 2] = het_probs[curr_hets] * curr_hets * (curr_hets - 1.0)
			/ (4.0 * (curr_homr + 1.0) * (curr_homc + 1.0));
		sum += het_probs[curr_hets - 2];

		// 2 fewer heterozygotes for next iteration -> add one rare, one common homozygote
		curr_homr++;
		curr_homc++;
	}

	curr_hets = mid;
	curr_homr = (rare_copies - mid) / 2;
	curr_homc = genotypes - curr_hets - curr_homr;
	for (curr_hets = mid; curr_hets <= rare_copies - 2; curr_hets += 2) {
		het_probs[curr_hets + 2] = het_probs[curr_hets] * 4.0 * curr_homr * curr_homc
			/((curr_hets + 2.0) * (curr_hets + 1.0));
		sum += het_probs[curr_hets + 2];

		// add 2 heterozygotes for next iteration -> subtract one rare, one common homozygote 
		curr_homr--;
		curr_homc--;
	}

	for (i = 0; i <= rare_copies; i++)
		het_probs[i] /= sum;

	//  p-value calculation for p_hwe 
	for (i = 0; i <= rare_copies; i++)
	{
		if (het_probs[i] > het_probs[obs_hets])
			continue;
		p_hwe += het_probs[i];
	}

	p_hwe = p_hwe > 1.0 ? 1.0 : p_hwe;
	free(het_probs);

	return p_hwe;
}

double normdist(double z) {
	double sqrt2pi =  2.50662827463;
	double t0, z1, p0;

	t0 = 1 / (1 +  0.2316419 *  fabs(z));
	z1 =  exp( -0.5 * z*z ) / sqrt2pi;
	p0 = z1 * t0
		* (0.31938153 +
		t0 * ( -0.356563782 +
		t0 * ( 1.781477937 +
		t0 * ( -1.821255978 +
		1.330274429 * t0))));

	return z >= 0 ? 1 - p0 : p0 ;
}

// Calculate correlation matrix for X
// Skip intercept
BOOL check_vif(MATRIX_T X, int x_rows, int x_cols) {
	int p, q, i, j, k;
	VECTOR_T vif_m = NULL;
	MATRIX_T vif_c = NULL, vif_c_inverse = NULL;

	p = x_rows;
	if (p < 2) return FALSE;

	q = x_cols - 1;
	if ( q < 2 ) return TRUE;

	// Assign the vector/matrix arrays required by this function.
	vif_m = create_vector(q, 0);
	vif_c = create_matrix(q, q, 0);
	vif_c_inverse = create_matrix(q, q, 0);
	if (vif_m == NULL || vif_c == NULL || vif_c_inverse == NULL) {
		if (vif_m != NULL) free(vif_m); 
		if (vif_c != NULL) free_matrix(&vif_c); 
		if (vif_c_inverse != NULL) free_matrix(&vif_c_inverse);
	}

	for (i = 0; i < p; i++)
		for (j = 0; j < q; j++)
			vif_m[j] += X[i][j + 1];

	for (j = 0; j < q; j++)
		vif_m[j] /=  p;

	for (i = 0; i < p; i++)
		for (j = 0; j < q; j++)
			for (k = j; k < q; k++)
				vif_c[j][k] += (X[i][j+1] - vif_m[j]) * (X[i][k+1] - vif_m[k]);

	for (j = 0; j < q; j++)
		for (k = j; k < q; k++) 
			vif_c[j][k] /= (p - 1);

	for (j = 0; j < q; j++) {
		for (k= j + 1; k < q; k++) {
			vif_c[j][k] /= sqrt(vif_c[j][j] * vif_c[k][k] );
			vif_c[k][j] = vif_c[j][k];

			// individual element > 0.999, skipping VIF
			if ( vif_c[k][j] > 0.999) {
				if (vif_m != NULL) free(vif_m); 
				if (vif_c != NULL) free_matrix(&vif_c);
				if (vif_c_inverse != NULL) free_matrix(&vif_c_inverse);
				return FALSE;
			}
		}		
	}

	for (j = 0; j < q; j++) {
		if (vif_c[j][j] == 0 || !realnum(vif_c[j][j])) {
			if (vif_m != NULL) free(vif_m); 
			if (vif_c != NULL) free_matrix(&vif_c);
			if (vif_c_inverse != NULL) free_matrix(&vif_c_inverse);
			return FALSE;
		}
		vif_c[j][j] = 1;
	}

	// Get inverse
	if (!svd_inverse(vif_c, vif_c_inverse, q, q)) return FALSE;

	// Calculate VIFs
	for (j = 0; j < q;j++) {
		// r^2 = 1 - 1/x where x is diagonal element of inverted
		// correlation matrix
		// As VIF = 1 / ( 1 - r^2 ) , implies VIF = x
		if (vif_c_inverse[j][j] > vif_threshold) {
			if (vif_m != NULL) free(vif_m); 
			if (vif_c != NULL) free_matrix(&vif_c);
			if (vif_c_inverse != NULL) free_matrix(&vif_c_inverse);
			return FALSE;		
		}
	}

	if (vif_m != NULL) free(vif_m); 
	if (vif_c != NULL) free_matrix(&vif_c);
	if (vif_c_inverse != NULL) free_matrix(&vif_c_inverse);

	return TRUE;
}

int svd_inverse(MATRIX_T u, MATRIX_T r, int nRows, int nCols) {
	int i, j, k, n = nRows;
	const double eps = DBL_EPSILON;
	double wmax = 0, wmin;
	VECTOR_T w = NULL;
	MATRIX_T v = NULL;

	// Check that matrix dimensions agree.
	if (nRows == 0 || nRows != nCols) return 0;

	// Allocate the local matrix & vector.
	w = create_vector(n, 0);
	v = create_matrix(n, n, 0);
	if (w == NULL || v == NULL) {
		if (w != NULL) free(w); w = NULL;
		if (v != NULL) free_matrix(&v); v = NULL;
	}

	if(!svdcmp(u, nRows, nCols, w, v)) return 0;

	// Look for singular values
	for (i = 0; i < n; i++) wmax = w[i] > wmax ? w[i] : wmax;
	wmin = wmax * eps;
	wmin = wmin < eps ? eps : wmin; 
	for (i = 0; i < n; i++) w[i] = w[i] < wmin ? 0 : 1/w[i];

	// row U * 1/w
	for (i = 0; i < n; i++) {
		for (j = 0; j < n; j++) u[i][j] = u[i][j] * w[j];
	}

	// [nxn].[t(v)]
	zero_matrix(r, n, n);
	for (i = 0; i < n; i++)
		for (j = 0; j < n; j++)
			for (k = 0; k < n; k++)
				r[i][j] += u[i][k] * v[j][k];

	if (w != NULL) free(w); w = NULL;
	if (v != NULL) free_matrix(&v); v = NULL;

	return 1;
}

int svdcmp(MATRIX_T a, int a_rows, int a_cols, VECTOR_T w, MATRIX_T v) {
	BOOL flag;
	int i, its ,j , jj, k, l = 0, nm = 0;
	double anorm, c, f, g, h, s, scale, x, y, z;
	double volatile temp;
	int m, n;
	VECTOR_T rv1 = NULL;

	m = a_rows;
	if (m == 0) return 0;
	n = a_cols;
	rv1 = create_vector(n, 0);
	if (rv1 == NULL) return 0;

	g=scale=anorm=0.0;
	for (i=0;i<n;i++) {
		l=i+2;
		rv1[i]=scale*g;
		g=s=scale=0.0;
		if (i < m) {
			for (k=i;k<m;k++) scale += fabs(a[k][i]);
			if (scale != 0.0) {
				for (k=i;k<m;k++) {
					a[k][i] /= scale;
					s += a[k][i]*a[k][i];
				}
				f=a[i][i];
				g = -SIGN(sqrt(s),f);
				h=f*g-s;
				a[i][i]=f-g;
				for (j=l-1;j<n;j++) {
					for (s=0.0,k=i;k<m;k++) s += a[k][i]*a[k][j];
					f=s/h;
					for (k=i;k<m;k++) a[k][j] += f*a[k][i];
				}
				for (k=i;k<m;k++) a[k][i] *= scale;
			}
		}
		w[i]=scale *g;
		g=s=scale=0.0;
		if (i+1 <= m && i+1 != n) {
			for (k=l-1;k<n;k++) scale += fabs(a[i][k]);
			if (scale != 0.0) {
				for (k=l-1;k<n;k++) {
					a[i][k] /= scale;
					s += a[i][k]*a[i][k];
				}
				f=a[i][l-1];
				g = -SIGN(sqrt(s),f);
				h=f*g-s;
				a[i][l-1]=f-g;
				for (k=l-1;k<n;k++) rv1[k]=a[i][k]/h;
				for (j=l-1;j<m;j++) {
					for (s=0.0,k=l-1;k<n;k++) s += a[j][k]*a[i][k];
					for (k=l-1;k<n;k++) a[j][k] += s*rv1[k];
				}
				for (k=l-1;k<n;k++) a[i][k] *= scale;
			}
		}
		anorm=MAX(anorm,(fabs(w[i])+fabs(rv1[i])));
	}
	for (i=n-1;i>=0;i--) {
		if (i < n-1) {
			if (g != 0.0) {
				for (j=l;j<n;j++)
					v[j][i]=(a[i][j]/a[i][l])/g;
				for (j=l;j<n;j++) {
					for (s=0.0,k=l;k<n;k++) s += a[i][k]*v[k][j];
					for (k=l;k<n;k++) v[k][j] += s*v[k][i];
				}
			}
			for (j=l;j<n;j++) v[i][j]=v[j][i]=0.0;
		}
		v[i][i]=1.0;
		g=rv1[i];
		l=i;
	}
	for (i=MIN(m,n)-1;i>=0;i--) {
		l=i+1;
		g=w[i];
		for (j=l;j<n;j++) a[i][j]=0.0;
		if (g != 0.0) {
			g=1.0/g;
			for (j=l;j<n;j++) {
				for (s=0.0,k=l;k<m;k++) s += a[k][i]*a[k][j];
				f=(s/a[i][i])*g;
				for (k=i;k<m;k++) a[k][j] += f*a[k][i];
			}
			for (j=i;j<m;j++) a[j][i] *= g;
		} else for (j=i;j<m;j++) a[j][i]=0.0;
		++a[i][i];
	}
	for (k=n-1;k>=0;k--) {
		for (its=0;its<30;its++) {
			flag=TRUE;
			for (l=k;l>=0;l--) {
				nm=l-1;
				temp=fabs(rv1[l])+anorm;
				if (temp == anorm) {
					flag=FALSE;
					break;
				}
				temp=fabs(w[nm])+anorm;
				if (temp == anorm) break;
			}
			if (flag) {
				c=0.0;
				s=1.0;
				for (i=l;i<k+1;i++) {
					f=s*rv1[i];
					rv1[i]=c*rv1[i];
					temp = fabs(f)+anorm;
					if (temp == anorm) break;
					g=w[i];
					h=pythag(f,g);
					w[i]=h;
					h=1.0/h;
					c=g*h;
					s = -f*h;
					for (j=0;j<m;j++) {
						y=a[j][nm];
						z=a[j][i];
						a[j][nm]=y*c+z*s;
						a[j][i]=z*c-y*s;
					}
				}
			}
			z=w[k];
			if (l == k) {
				if (z < 0.0) {
					w[k] = -z;
					for (j=0;j<n;j++) v[j][k] = -v[j][k];
				}
				break;
			}
			if (its == 30) {
				error("no convergence in 30 svdcmp iterations");
				free(rv1);
				return 0;
			}
			x=w[l];
			nm=k-1;
			y=w[nm];
			g=rv1[nm];
			h=rv1[k];
			f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
			g=pythag(f,1.0);
			f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
			c=s=1.0;
			for (j=l;j<=nm;j++) {
				i=j+1;
				g=rv1[i];
				y=w[i];
				h=s*g;
				g=c*g;
				z=pythag(f,h);
				rv1[j]=z;
				c=f/z;
				s=h/z;
				f=x*c+g*s;
				g=g*c-x*s;
				h=y*s;
				y *= c;
				for (jj=0;jj<n;jj++) {
					x=v[jj][j];
					z=v[jj][i];
					v[jj][j]=x*c+z*s;
					v[jj][i]=z*c-x*s;
				}
				z=pythag(f,h);
				w[j]=z;
				if (z) {
					z=1.0/z;
					c=f*z;
					s=h*z;
				}
				f=c*g+s*y;
				x=c*y-s*g;
				for (jj=0;jj<m;jj++) {
					y=a[jj][j];
					z=a[jj][i];
					a[jj][j]=y*c+z*s;
					a[jj][i]=z*c-y*s;
				}
			}
			rv1[l]=0.0;
			rv1[k]=f;
			w[k]=x;
		}
	}
	if (rv1 != NULL) free(rv1);
	
	return 1;
}

int mult_matrix(MATRIX_T a, int ar, int ac, MATRIX_T b, int br, int bc, MATRIX_T c) {
	int i, j, k;

	// Check the matrix
	if (ar == 0 || br == 0 || ac != br) 
		return 0;

	for (i = 0; i <ar; i++)	
		for (j = 0; j < bc; j++) 
			for (k = 0; k < ac; k++) 
				c[i][j] += a[i][k] * b[k][j];

	return 1;
}

enum MATHS_RETURN svbksb(MATRIX_T u, int u_rows, int u_cols, VECTOR_T w, MATRIX_T v, VECTOR_T b, VECTOR_T x) {
	int jj,j,i;
	double s;
	int m = u_rows, n = u_cols;
	VECTOR_T tmp = create_vector(n, 0);

	if (tmp == NULL) return MEMORY_FAILURE;

	for (j = 0; j < n; j++) {
		s=0.0;
		if (w[j] != 0.0) {
			for (i = 0; i < m; i++) s += u[i][j] * b[i];
			s /= w[j];
		}
		tmp[j]=s;
	}

	for (j = 0; j < n; j++) {
		s=0.0;
		for (jj = 0; jj < n; jj++) s += v[j][jj] * tmp[jj]; 
		x[j] = s;
	}
	if (tmp != NULL) free(tmp);

	return MATHS_SUCCESS;
}

long getPairwiseTests(long n) { 
	if (n == 0) return 0;
	else return (long)((n*(n-1))/2); 
}


long getPairwiseTestsPerNode(long nElements, long nNodes) {
	return (long) (nElements * (nElements - 1)) / (2 *nNodes); 
}

// Algorithm by David Hoyle (david.hoyle@manchester.ac.uk), 16/08/2008
BOOL getDHLoopLimits(int N, int nodes, int mpi_rank, int *indices) {
	unsigned long nTotalTests, nTestsPerNode;
	int i1 = 1, i2, k;

	indices[0] = indices[1] = -1;
	nTotalTests = getPairwiseTests(N);
	nTestsPerNode = getPairwiseTestsPerNode(N, nodes);

	for (k = 1; k <= nodes && (k - 1) <= mpi_rank; k++) {
		if (k < nodes) {
			// define quadratic coefficients
			double b = -(((double) N) - 0.5);
			double c = nTestsPerNode + (N * (i1 - 1)) - (0.5 * i1 * (i1 - 1));
			double discriminant = (b * b) - (2.0 * c);
			
			if (discriminant >= 0.0) {
				double temp = -(b) - sqrt(discriminant);

				// Assume negative branch.
				i2 = (int) temp;
			} else {
				// Complex solution
				return FALSE;
			} 
		} else i2 = N;

		// Record the loop limits.
		indices[0] = i1 - 1;
		indices[1] = i2;

		i1 = i2 + 1;
	}

	return TRUE;
}

double getChisqPvalue(double chisq) { return normdist(-(sqrt(chisq))) * 2; }

// Update the permutation counts for 1-D array of results based on the PLINK schema.
int perm_update(int size, double *results, double *originals, int *R, int *maxR) {
	double mx = 0;
	int i;

	if (size <= 0 || results == NULL || originals == NULL || R == NULL || maxR == NULL) 
		return 0;

	// Update number of successes
	for (i = 0; i < size; i++) if (results[i] >= originals[i] || !realnum(originals[i])) R[i]++;

	// Find data-set maximum.
	for (i = 0; i < size; i++)
		if (realnum(results[i]))
			if (results[i] > mx) mx = results[i];

	// Max(T) permutation -- compare against best
	for (i = 0; i < size; i++) if (mx >= originals[i] || !realnum(originals[i])) maxR[i]++;

	return 1;
}

// Based on the Fisher-Yates algorithm
void shuffle(int *src, int nSrc) {
	int n;

	if (nSrc <= 0 && src != NULL) return;


	n = nSrc;
	while (n > 1)  {
		int temp;
		int k = get_rand(n);// 0 <= k < n.

		n--;			// n is now the last pertinent index;
		temp = src[n];	// swap array[n] with array[k] (does nothing if k == n).
		src[n] = src[k];
		src[k] = temp;
	}
}

unsigned long get_tasks_per_processor(unsigned int numprocs, unsigned long nSample) {
	double average = (double) nSample /(double) numprocs;
	return (unsigned long) floor(average);
}

int _get_tasks_per_processor(int numprocs, int nSample) {
	double average = (double) nSample /(double) numprocs;
	return (int) ceil(average);
}

// Method provided by Wolf, see http://slack.ser.man.ac.uk/theory/association_odds.html
// Based on asymptotic (large sample) approximation.
int getOddsRatio(int A1, int A2, int U1, int U2, double *OR, double *lowerCI, double *upperCI, double zt) {
	if (OR == NULL || lowerCI == NULL || upperCI == NULL) return 0;

	*OR = *lowerCI = *upperCI = 0;
	if (A1 != 0 && A2 != 0 && U1 != 0 && U2 != 0) {
		double one = 1, log_odds;
		double sum_of_reciprocals = (one /(double) A1) + (one /(double) A2) + (one /(double) U1) + (one /(double) U2);
		double limit = zt * sqrt(sum_of_reciprocals);

		*OR = (double) (A1 * U2) / (double) (U1 * A2);
		log_odds = log(*OR);
		*lowerCI = exp(log_odds - limit);
		*upperCI = exp(log_odds + limit);
	} else 
		return 0;


	return 1;
}